


#include "sl_newton.h"


// 这里实现牛顿下降法。   主要是实现 Broyden算法。《统计学习方法》李航 第二版 P446。



/// <summary>
/// x：自变量，返回目标函数在x处的值。
/// </summary>
typedef float (*target_function)(vec x);

/// <summary>
///  x：自变量，gf：函数在x处的偏导数
/// </summary>
typedef void (*function_grad)(vec x, vec gf);


/// <summary>
/// 根据三组 x，y的坐标，用二次曲线拟合，估计极小值点。 parabola:抛物线。
/// </summary>
/// <param name="x0"></param>
/// <param name="y0"></param>
/// <param name="x1"></param>
/// <param name="y1"></param>
/// <param name="x2"></param>
/// <param name="y2"></param>
/// <returns>估计的极小值的x坐标</returns>
float get_parabola_min(float x1, float y1, float x2, float y2, float x3, float y3)
{
	// 用二次曲线拟合计算lamda的极小值点。
	float a, b ;// , c;
	/*
//（x1,y1),(x2,y2),(x3,y3)可以确定一个二次函数：
//y = y3*(x-x1)(x-x2)/[(x3-x1)(x3-x2)]+y2*(x-x1)(x-x3)/[(x2-x1)(x2-x3)]+y1*(x-x2)(x-x3)/[(x1-x2)(x1-x3)]
//  (已知x1=0,x2=0.1,x3=0.2) 得：
a =   y3/[(x3-x1)(x3-x2)]        +y2/[(x2-x1)(x2-x3)]        +y1/[(x1-x2)(x1-x3)]
b =-{ y3*(x1+x2)/[(x3-x1)(x3-x2)]+y2*(x1+x3)/[(x2-x1)(x2-x3)]+y1*(x2+x3)/[(x1-x2)(x1-x3)] }
c =   y3*(x1*x2)/[(x3-x1)(x3-x2)]+y2*(x1*x3)/[(x2-x1)(x2-x3)]+y1*(x2*x3)/[(x1-x2)(x1-x3)]

(x3-x1)(x3-x2) = 0.02,  (x2-x1)(x2-x3) = -0.01,  (x1-x2)(x1-x3) = 0.02,
(x1+x2) = 0.1    (x1+x3) = 0.2   (x2+x3) = 0.3
(x1*x2) = 0      (x1*x3) = 0     (x2*x3) = 0.02
a = (y3 - 2*y2 + y1)/0.02
b = (y3 - 4*y2 + 3*y1)/0.2
c = y1               //  下面对 abc进行分类讨论，计算x的极小值点。注意：x>0.
		a = (y3 - 2 * y2 + y1) / 0.02;
		b = (y3 - 4 * y2 + 3 * y1) / 0.2;
		c = y1;
		*/

	float x3_x1_x3_x2 = ((x3 - x1) * (x3 - x2));
	float x2_x1_x2_x3 = ((x2 - x1) * (x2 - x3));
	float x1_x2_x1_x3 = ((x1 - x2) * (x1 - x3));
	a =  y3 / x3_x1_x3_x2              + y2 / x2_x1_x2_x3             + y1 / x1_x2_x1_x3;
	b = -(y3 * (x1 + x2) / x3_x1_x3_x2 + y2 * (x1 + x3) / x2_x1_x2_x3 + y1 * (x2 + x3) / x1_x2_x1_x3);
//	c =  y3 * (x1 * x2) / x3_x1_x3_x2  + y2 * (x1 * x3) / x2_x1_x2_x3 + y1 * (x2 * x3) / x1_x2_x1_x3;
	float lamdak = 1;
	if (ABS(a) < Epsilon) {  //这里拟合出来是线性函数。
		if (b > 0)   // 线性上升……  只能期待后面有拐点了。 一般应该不会遇到这种情况
			lamdak = x1 + (x3 - x1) * 10;
		else // 线性下降。
			lamdak = x1 + (x3 - x1) * 3;
	}
	else
	{
		float x_ex = -b / (2 * a); // 极值点。
		if (a > 0) {  // 开口向上，有极小值点
			if (x_ex > x1)  // 极小值点在右边，计算该点的函数值，如果比其他点都小，就选它，否则选三个点中最小的。
				lamdak = x_ex;
			else // (x_ex < x1)   // 极小值点在左边，说明x1,x2都是在增加。说明步长可能太大，最小值可能在(x0,x1)之间。
				lamdak = (x1 + x2) / 2;
		}
		else  // 开口朝下，有极大值点。
		{
			if (x_ex > x1)  // 极值点在右边，说明x1,x2都是在增加。如果要比现在的点小，按照拟合函数，就得大步跨过极值点：
				lamdak = x_ex + (x_ex - x1) + (x3 - x1) * 3;
			else //  (x_ex < x1)   // 极值点在左边，说明后面在加速减小。
				lamdak = x1 + (x3 - x1) * 3;
		}
	}
	//

	return lamdak;
}



/// <summary>
/// 这里实现Sherman_Morrison 公式：列矢u,v,的转置是uT,vT,   A*iA = I >>  (A+u*vT)^(-1) = iA - (iA*u*vT*iA)/(1+vT*iA*u) 这个公式求逆，计算量比较小
/// "A"  已知正定矩阵，其实A不需要参与运算……，因此没有这个参数
/// </summary>
/// <param name="iA">已知A的逆</param>
/// <param name="u">已知列矢</param>
/// <param name="v">已知列矢</param>
/// <param name="iAu">中间的计算变量，无需初始化</param>
/// <param name="vtiA">中间的计算变量，无需初始化</param>
/// <param name="iAuvt">待求矩阵(A+u*vT)的逆</param>
/// <returns>是否运算成功。只要矩阵元素的尺寸匹配都可以成功运行。</returns>
bool mat_Sherman_Morrison( mat iA, vec u, vec v, vec iAu, vec vtiA, mat iAuvt)
{
	mat_mul_vec(iA, u, iAu);
	matT_mul_vec(iA, v, vtiA);

	float a = 1 + vec_dotProduct(v, iAu);
	mat_vec_mul_vecT(iAuvt, iAu, vtiA);
	mat_mul_real(iAuvt, 1.0f / a, iAuvt);
	mat_sub(iA, iAuvt, iAuvt);
	return true;
}
//

/// <summary>
/// 这里实现 对称矩阵的Sherman_Morrison 公式，对于对称的情况可以减少计算量：列矢u,的转置是uT,   A*iA = I >>  (A+u*uT)^(-1) = iA - (iA*u*uT*iA)/(1+uT*iA*u) 
/// "A"  已知正定矩阵，其实A不需要参与运算……，因此没有这个参数。
/// 这个公式求逆，计算量比较小
/// </summary>
/// <param name="iA">已知A的逆</param>
/// <param name="u">已知列矢</param>
/// <param name="iAu">中间的计算变量，无需初始化</param>
/// <param name="iAuvt">待求矩阵(A+u*vT)的逆</param>
/// <returns>是否运算成功。只要矩阵元素的尺寸匹配都可以成功运行。</returns>
bool mat_Symmetric_Sherman_Morrison_add(mat iA, vec u, vec iAu,  mat iAuvt)
{
	mat_mul_vec(iA, u, iAu);
	float a = 1.0f / (1 + vec_dotProduct(u, iAu));
	real_t** iA_data = iA.data;
	real_t** iAuvt_data = iAuvt.data;
	real_t* iAu_data = iAu.data;
	for (int i = 0; i < iA.rows; i++) {
		for (int j = i; j < iA.cols; j++) {
			iAuvt_data[i][j] = iA_data[i][j] - iAu_data[i] * iAu_data[j] * a;
			iAuvt_data[j][i] = iAuvt_data[i][j];
		}
	}
	return true;
}
//

/// <summary>
/// 这里实现 对称矩阵的Sherman_Morrison 公式，对于对称的情况可以减少计算量：列矢u,的转置是uT,   A*iA = I >>  (A-u*uT)^(-1) = iA + (iA*u*uT*iA)/(1-uT*iA*u) 
/// "A"  已知正定矩阵，其实A不需要参与运算……，因此没有这个参数
/// </summary>
/// <param name="iA">已知A的逆</param>
/// <param name="u">已知列矢</param>
/// <param name="iAu">中间的计算变量，无需初始化</param>
/// <param name="iAuvt">待求矩阵(A+u*vT)的逆</param>
/// <returns>是否运算成功。只要矩阵元素的尺寸匹配都可以成功运行。</returns>
bool mat_Symmetric_Sherman_Morrison_sub(mat iA, vec u, vec iAu, mat iAuvt)
{
	mat_mul_vec(iA, u, iAu);
	float a = 1.0f / (1 - vec_dotProduct(u, iAu));
	real_t** iA_data = iA.data;
	real_t** iAuvt_data = iAuvt.data;
	real_t* iAu_data = iAu.data;
	for (int i = 0; i < iA.rows; i++) {
		for (int j = i; j < iA.cols; j++) {
			iAuvt_data[i][j] = iA_data[i][j] + iAu_data[i] * iAu_data[j] * a;
			iAuvt_data[j][i] = iAuvt_data[i][j];
		}
	}
	return true;
}
//



#include<math.h>

/// <summary>
/// 拟牛顿法，计算目标函数f的极小点xmin,返回迭代次数。
/// </summary>
/// <param name="f">目标函数</param>
/// <param name="fg">目标函数的导数函数</param>
/// <param name="err">精度要求</param>
/// <param name="x_min">极小值点，同时提供初始点</param>
int BFGS(target_function f, function_grad fg, float err, vec xmin)
{
	int size = xmin.size;  // 自变量的元素个数
	int k = 0;
	vec xk, xk1, gk, gk1, pk, vt2, vt3, vt4;
	mat bt, ibt, ibk, ibk1, bk, bk1;
	{
		vec_allocate(&xk, size);
		vec_allocate(&xk1, size);
		vec_allocate(&gk, size);
		vec_allocate(&gk1, size);
		vec_allocate(&pk, size);
		vec_allocate(&vt2, size);
		vec_allocate(&vt3, size);
		vec_allocate(&vt4, size);

		mat_allocate(&bt, size, size);
		mat_allocate(&ibt, size, size);
		mat_allocate(&bk, size, size);
		mat_allocate(&bk1, size, size);
		mat_allocate(&ibk, size, size);
		mat_allocate(&ibk1, size, size);
	}

	// 初始化矩阵Bk为单位矩阵,将初始点赋值给xk。
	mat_eye(bk);
	mat_eye(ibk);
	vec_copyfrom(xk, xmin);
	fg(xk, gk);
	float gk_length = vec_length(gk);
	if (gk_length < err) {
		vec_copyto(xk, xmin);
		return k;
	}

	// 下面按照书上的算法B.3 的流程走。

	float y1 = f(xk), y2, y3;
	float x1 = 0, x2 = 0.1f, x3 = 0.2f;   // x1 < x2 < x3 ,要求的极小值点在x1右边，(xmin > x1)


	while (k < 50)  // 迭代次数。
	{
		//printfs("  \r\n  %d \r\n   ", k);
		//vec_print(xk);
		//vec_print(gk);

		k++;
		mat_mul_vec(ibk, gk, pk);
		for (int i = 0; i < size; i++) {
			vt2.data[i] = xk.data[i] - x2 * pk.data[i];
			vt3.data[i] = xk.data[i] - x3 * pk.data[i];
		}
		y2 = f(vt2), y3 = f(vt3);

		// 利用三个坐标点，二次曲线拟合，估计极小值点。
		float x4 = get_parabola_min(x1, y1, x2, y2, x3, y3);
		// 计算 lamdak 处的值，然后在 x2,x3,x4 中选函数值最小的点，作为xk1。
		for (int i = 0; i < size; i++) {
			vt4.data[i] = xk.data[i] - x4 * pk.data[i];
		}
		float y4 = f(vt4);
		// 计算xk1,刷新y1
		if (y2 <= y3 && y2 <= y4) {
			vec_copyfrom(xk1, vt2);
			y1 = y2;
		}
		if (y3 <= y2 && y3 <= y4) {
			vec_copyfrom(xk1, vt3);
			y1 = y3;
		}
		if (y4 <= y2 && y4 <= y3) {
			vec_copyfrom(xk1, vt4);
			y1 = y4;
		}

		// 计算gk1
		fg(xk1, gk1);
		gk_length = vec_length(gk1);
		if (gk_length < err) {
			vec_copyto(xk, xmin);
			return k;
		}

		//刷新bk,ibk;  得先得到 vt2 = dk = xk1-xk,  vt3 = yk = gk1 - gk; 
		vec dk = vt2;
		vec yk = vt3;
		vec bkdk = vt4;
		vec_sub(xk1, xk, dk);
		vec_sub(gk1, gk, yk);
		mat_mul_vec(bk, dk, bkdk); //

		float yd = vec_dotProduct(yk, dk);     // 这里有可能出现负号的情况。 要做相应的处理
		float dbd = vec_dotProduct(dk, bkdk);
		int sign_yd = yd > 0 ? 1 : -1;
		int sign_dbd = dbd > 0 ? 1 : -1;
		yd *= sign_yd;
		dbd *= sign_dbd;

		// printfs(" \r\n  yd = %f    dbd = %f \r\n  ", yd, dbd);

		vec_mul(yk, 1.0f / sqrtf(yd), yk);
		vec_mul(bkdk, 1.0f / sqrtf(dbd), bkdk);

		// 更新bk1
		float** bk1_data = bk1.data;
		float** bk_data = bk.data;
		float** bt_data = bt.data;
		float* yk_data = yk.data;
		float* bkdk_data = bkdk.data;
		for (int i = 0; i < size; i++) {
			for (int j = i; j < size; j++) {
				bt_data[i][j] = bk_data[i][j] + yk_data[i] * yk_data[j] * sign_yd;
				bt_data[j][i] = bt_data[i][j];

				bk1_data[i][j] = bt_data[i][j] - bkdk_data[i] * bkdk_data[j] * sign_dbd;
				bk1_data[j][i] = bk1_data[i][j];
			}
		}
		if (sign_yd > 0) {
			mat_Symmetric_Sherman_Morrison_add(ibk, yk, vt2, ibt);  // 现在bt 是 ibt 的逆。
		}
		else
		{
			mat_Symmetric_Sherman_Morrison_sub(ibk, yk, vt2, ibt);  // 现在bt 是 ibt 的逆。
		}
		if (sign_dbd > 0) {
			mat_Symmetric_Sherman_Morrison_sub(ibt, bkdk, vt2, ibk1);
		}
		else
		{
			mat_Symmetric_Sherman_Morrison_add(ibt, bkdk, vt2, ibk1);
		}

		//mat_mul(bk1, ibk1, bk);
		//mat_print(bk1);
		//mat_print(ibk1);
		//mat_print(bk);
		mat_copyfrom(bk, bk1);
		mat_copyfrom(ibk, ibk1);
		vec_copyfrom(gk, gk1);
		vec_copyfrom(xk, xk1);
	}

	vec_copyto(xk, xmin);
	{
		vec_deallocate(&xk);
		vec_deallocate(&xk1);
		vec_deallocate(&gk);
		vec_deallocate(&gk1);
		vec_deallocate(&pk);
		vec_deallocate(&vt2);
		vec_deallocate(&vt3);
		vec_deallocate(&vt4);
		mat_deallocate(&bt);
		mat_deallocate(&ibt);
		mat_deallocate(&ibk);
		mat_deallocate(&ibk1);
		mat_deallocate(&bk);
		mat_deallocate(&bk1);
	}
	return k;

}






float f1(float x)
{
	float a = 0.4, b = -0.9, c = 4;
	return a * x * x + b * x + c;
}

void testmat_Sherman_Morrison(void)
{
	int row = 13;
	int col = 13;
	int Length = 13;
	mat m;
	mat n;
	mat im;
	mat in;
	mat uvt;
	mat eye;
	vec A;
	vec B;
	vec C;
	vec D;
	vec E;
	vec_allocate(&A, Length);
	vec_allocate(&B, Length);
	vec_allocate(&C, Length);
	vec_allocate(&D, Length);
	vec_allocate(&E, Length);

	mat_allocate(&m, row, col);
	mat_allocate(&n, row, row);
	mat_allocate(&im, row, col);
	mat_allocate(&in, row, col);
	mat_allocate(&uvt, row, row);
	mat_allocate(&eye, row, row);

	vec_randomize(A, 314159, 0, 1);
	vec_randomize(B, 31415, 0, 1);
	mat_randomize(m, 314159, -2, 2);
	//for (int i = 0; i < row; i += 1) {
	//	for (int j = 0; j < col; j += 1) {
	//		m.data[j][i] = m.data[i][j];
	//	}
	//}

	mat_vec_mul_vecT(uvt, A, B);
	mat_add(m, uvt, n);
	mat_print(uvt);
	mat_inversion(m, im, Epsilon);
	mat_print(im);

	mat_Sherman_Morrison(im, A, B, C, D, in);
	vec_print(A);
	vec_print(B);
	mat_mul(n, in,eye);
	mat_print(n);
	mat_print(in);
	mat_print(eye);

	mat_deallocate(&m);  //不要求m矩阵是方阵。A尺寸与m矩阵一样。H是和行数一样的方阵。
	mat_deallocate(&n);
	mat_deallocate(&im);
	mat_deallocate(&in);
	mat_deallocate(&uvt);
	mat_deallocate(&eye);
	vec_deallocate(&A);
	vec_deallocate(&B);
	vec_deallocate(&C);
	vec_deallocate(&D);
	vec_deallocate(&E);

}

void testmat_Symmetric_Sherman_Morrison_add(void)
{
	int row = 13;
	int col = 13;
	int Length = 13;
	mat m;
	mat n;
	mat im;
	mat in;
	mat uvt;
	mat eye;
	vec A;
	vec B;
	vec C;
	vec D;
	vec E;
	vec_allocate(&A, Length);
	vec_allocate(&B, Length);
	vec_allocate(&C, Length);
	vec_allocate(&D, Length);
	vec_allocate(&E, Length);

	mat_allocate(&m, row, col);
	mat_allocate(&n, row, row);
	mat_allocate(&im, row, col);
	mat_allocate(&in, row, col);
	mat_allocate(&uvt, row, row);
	mat_allocate(&eye, row, row);

	vec_randomize(A, 314159, 0, 1);
//	vec_randomize(B, 31415, 0, 1);
	mat_randomize(m, 314159, -2, 2);
	for (int i = 0; i < row; i += 1) {
		for (int j = 0; j < col; j += 1) {
			m.data[j][i] = m.data[i][j];
		}
	}

	vec_copyto(A, B);
	mat_vec_mul_vecT(uvt, A, B);
	mat_add(m, uvt, n);
	mat_print(uvt);
	mat_inversion(m, im, Epsilon);
	mat_print(im);
	mat_Symmetric_Sherman_Morrison_add(im, A, C, in);


	//mat_Sherman_Morrison(im, A, B, C, D, in);

	vec_print(A);
	vec_print(B);
	mat_mul(n, in, eye);
	mat_print(n);
	mat_print(in);
	mat_print(eye);

	mat_deallocate(&m);  //不要求m矩阵是方阵。A尺寸与m矩阵一样。H是和行数一样的方阵。
	mat_deallocate(&n);
	mat_deallocate(&im);
	mat_deallocate(&in);
	mat_deallocate(&uvt);
	mat_deallocate(&eye);
	vec_deallocate(&A);
	vec_deallocate(&B);
	vec_deallocate(&C);
	vec_deallocate(&D);
	vec_deallocate(&E);

}

void testmat_Symmetric_Sherman_Morrison_sub(void)
{
	int row = 13;
	int col = 13;
	int Length = 13;
	mat m;
	mat n;
	mat im;
	mat in;
	mat uvt;
	mat eye;
	vec A;
	vec B;
	vec C;
	vec D;
	vec E;
	vec_allocate(&A, Length);
	vec_allocate(&B, Length);
	vec_allocate(&C, Length);
	vec_allocate(&D, Length);
	vec_allocate(&E, Length);

	mat_allocate(&m, row, col);
	mat_allocate(&n, row, row);
	mat_allocate(&im, row, col);
	mat_allocate(&in, row, col);
	mat_allocate(&uvt, row, row);
	mat_allocate(&eye, row, row);

	vec_randomize(A, 314159, 0, 1);
	//	vec_randomize(B, 31415, 0, 1);
	mat_randomize(m, 314159, -2, 2);
	for (int i = 0; i < row; i += 1) {
		for (int j = 0; j < col; j += 1) {
			m.data[j][i] = m.data[i][j];
		}
	}

	vec_copyto(A, B);
	mat_vec_mul_vecT(uvt, A, B);
	mat_sub(m, uvt, n);
	mat_print(uvt);
	mat_inversion(m, im, Epsilon);
	mat_print(im);
	mat_Symmetric_Sherman_Morrison_sub(im, A, C, in);


	//mat_Sherman_Morrison(im, A, B, C, D, in);

	vec_print(A);
	vec_print(B);
	mat_mul(n, in, eye);
	mat_print(n);
	mat_print(in);
	mat_print(eye);

	mat_deallocate(&m);  //不要求m矩阵是方阵。A尺寸与m矩阵一样。H是和行数一样的方阵。
	mat_deallocate(&n);
	mat_deallocate(&im);
	mat_deallocate(&in);
	mat_deallocate(&uvt);
	mat_deallocate(&eye);
	vec_deallocate(&A);
	vec_deallocate(&B);
	vec_deallocate(&C);
	vec_deallocate(&D);
	vec_deallocate(&E);

}




mat a;  
vec y; // 已知a，y，求最小的x，使得a*x = y.
vec x_min; // 已知a，y，求最小的x，使得a*x = y.
vec x; // 已知a，y，求最小的x，使得a*x = y.

float calcul_co(vec x) {
	vec b;
	vec err;
	vec_allocate(&b,y.size);
	vec_allocate(&err,y.size);

	mat_mul_vec(a, x, b);
	vec_sub(y, b, err);
	float s = vec_length_square(err);
	vec_deallocate(&b);
	vec_deallocate(&err);
	return s;
}

void grad_co(vec x, vec gx) {
	//mat a;
	//vec y; // 已知a，y，求最小的x，使得a*x = y.

	// 矩阵 M ，列矢v  d(M*v) = dM * v + M * dv >> 
	// 计算 d (y - a * x)^2  =  2 * (y - a * x) * d(y - a * x)   >>  d(y - a * x) 是列矢 = dy - ( da*x + a*dx ) >> 
	//   d (y - a * x)^2  =  2 * (y - a * x) * (-a*dx)    就是 矢量(y - a * x) 和矢量(a*dx)的点乘。 导数就很显然了：-2 * a^T * (y - a * x). 
	vec b;
	vec err;
	vec_allocate(&b, y.size);
	vec_allocate(&err, y.size);

	mat_mul_vec(a, x, b);
	vec_sub(y, b, err);
	matT_mul_vec(a, err, gx);
	vec_mul(gx, -2, gx);
	vec_deallocate(&b);
	vec_deallocate(&err);
}


void init_data()
{
	//m*x=v;已知m与v，求x。
	int row = 11;
	int col = 5;
	mat_allocate(&a, row, col);
	vec_allocate(&y, row);
	vec_allocate(&x, col);
	vec_allocate(&x_min, col);

	mat_randomize(a, 314159, -2, 2);
	vec_randomize(y, 27182, -2, 2);

	mat_LinearEquation_LinearLeastSquare(a, x_min, y, true, Epsilon);
	
	vec_init(x, 100);

	mat_print(a);
	vec_print(y);
	vec_print(x_min);
	vec_print(x);

}

void deinit_data()
{
	mat_deallocate(&a);
	vec_deallocate(&y);
	vec_deallocate(&x);
	vec_deallocate(&x_min);
}

void test_grad() {
	init_data();
//	int row = 11;
	int col = 5;
	vec v0,v,r;
	vec_allocate(&v0, col);
	vec_allocate(&r, col);
	vec_allocate(&v, col);

	vec_randomize(v0, 28, -2, 2);
	vec_randomize(r, 225, -2, 2);

	float a  = calcul_co(v0);
	float da;
	for (int i = 1; i < 4; i++) {
		 

		vec_mul(r, 0.01 * i, v);
		vec_add(v, v0, v);
		da = calcul_co(v) - a;
		a = calcul_co(v);

		grad_co(v, x);
		float g = vec_dotProduct(x, r) * 0.01;

		printfs(" \r\n  %d,  a = %f ,da = %f   g = %f ", i, a, da, g);
	}


	vec_deallocate(&v0);
	vec_deallocate(&v);
	vec_deallocate(&r);


	deinit_data();
}

void test_bgfs()
{
	init_data();

	int t = BFGS(calcul_co, grad_co, 0.000001f, x);

	vec_print(x);
	printfs("  \r\n  BFGS ==>> %d \r\n  ", t);
	deinit_data();

}


void test_bfgs_sub() 
{
	// 测试二次函数求极值
	{
		float x0, y0, x1, y1, x2, y2;
		x0 = 1;
		x1 = 2;
		x2 = 3;
		y0 = f1(x0);
		y1 = f1(x1);
		y2 = f1(x2);
		float x = get_parabola_min(x0, y0, x1, y1, x2, y2);
		printfs(" \r\n  (%f  , %f)  (%f  , %f)  (%f  , %f)  ", x0, y0, x1, y1, x2, y2);
		printfs(" \r\n x_min =  %f  y >> %f ", x, f1(x));
	}

	if(0)
	{
		printfs(" \r\n test  mat_Sherman_Morrison ");		
		testmat_Sherman_Morrison();

		printfs(" \r\n test  testmat_Symmetric_Sherman_Morrison_add  ");
		testmat_Symmetric_Sherman_Morrison_add();

		printfs(" \r\n test  testmat_Symmetric_Sherman_Morrison_sub  ");
		testmat_Symmetric_Sherman_Morrison_sub();

	}

//	test_grad();
	test_bgfs();
}


